Antonio Bartolomé Redondo
Daniel Buendía Ureña
# Importar las librerías necesarias
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt
import plotly.express as px
import sklearn
from sklearn.model_selection import train_test_split
from imblearn.over_sampling import SMOTE
from collections import Counter
from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import StandardScaler, OneHotEncoder
import category_encoders as ce
from sklearn.ensemble import RandomForestClassifier
from xgboost import XGBClassifier
from lightgbm import LGBMClassifier
from sklearn.metrics import confusion_matrix, accuracy_score, recall_score, precision_score, f1_score,\
roc_curve, precision_recall_curve
from sklearn.feature_selection import SelectFromModel
from sklearn import linear_model
from sklearn.linear_model import Lasso, LogisticRegression
from sklearn.metrics import plot_confusion_matrix
from sklearn.model_selection import GridSearchCV, RandomizedSearchCV
import scikitplot as skplt
from sklearn.inspection import PartialDependenceDisplay
import shap
import warnings
warnings.filterwarnings("ignore")
# Definición de funciones
def evaluate_model(ytest, ypred):
'''
La función permite calcular y mostrar las métricas básicas de la perfomance de un modelo de manera conjunta,
con el objetivo de que éste pueda ser fácilmente evaluado
'''
print("Accuracy:",accuracy_score(ytest, ypred))
print("Recall:",recall_score(ytest, ypred))
print("Precision:",precision_score(ytest, ypred))
print('F-Score: %.5f' % f1_score(ytest, ypred))
print('Confusion matrix: \n{}\n'.format(confusion_matrix(ytest, ypred)))
def get_feature_names(column_transformer):
"""Get feature names from all transformers.
Returns
-------
feature_names : list of strings
Names of the features produced by transform.
"""
# Remove the internal helper function
#check_is_fitted(column_transformer)
# Turn loopkup into function for better handling with pipeline later
def get_names(trans):
# >> Original get_feature_names() method
if trans == 'drop' or (
hasattr(column, '__len__') and not len(column)):
return []
if trans == 'passthrough':
if hasattr(column_transformer, '_df_columns'):
if ((not isinstance(column, slice))
and all(isinstance(col, str) for col in column)):
return column
else:
return column_transformer._df_columns[column]
else:
indices = np.arange(column_transformer._n_features)
return ['x%d' % i for i in indices[column]]
if not hasattr(trans, 'get_feature_names'):
# >>> Change: Return input column names if no method avaiable
# Turn error into a warning
warnings.warn("Transformer %s (type %s) does not "
"provide get_feature_names. "
"Will return input column names if available"
% (str(name), type(trans).__name__))
# For transformers without a get_features_names method, use the input
# names to the column transformer
if column is None:
return []
else:
return [name + "__" + f for f in column]
return [name + "__" + f for f in trans.get_feature_names()]
### Start of processing
feature_names = []
# Allow transformers to be pipelines. Pipeline steps are named differently, so preprocessing is needed
if type(column_transformer) == sklearn.pipeline.Pipeline:
l_transformers = [(name, trans, None, None) for step, name, trans in column_transformer._iter()]
else:
# For column transformers, follow the original method
l_transformers = list(column_transformer._iter(fitted=True))
for name, trans, column, _ in l_transformers:
if type(trans) == sklearn.pipeline.Pipeline:
# Recursive call on pipeline
_names = get_feature_names(trans)
# if pipeline has no transformer that returns names
if len(_names)==0:
_names = [name + "__" + f for f in column]
feature_names.extend(_names)
else:
feature_names.extend(get_names(trans))
return feature_names
# Lectura de los datos
data = pd.read_csv('data/pd_data_grouped_by_accident.csv', low_memory=False) # dtype=dtypes,
data
| C_MNTH | C_WDAY | C_HOUR | C_SEV | C_VEHS | C_CONF | C_RCFG | C_WTHR | C_RSUR | C_RALN | ... | V_TYPE_5 | V_TYPE_6 | V_TYPE_7 | V_TYPE_8 | V_TYPE_9 | P_SEX_F | P_SEX_M | P_SAFE | P_AGE | A_VAGE | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | 1 | 0 | 0 | 1 | 1 | 1 | 1 | 4 | 1 | ... | 0 | 1 | 0 | 0 | 0 | 1 | 0 | 1.0 | 17 | 8 |
| 1 | 1 | 1 | 0 | 0 | 1 | 1 | 1 | 1 | 5 | 1 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0.0 | 30 | 18 |
| 2 | 1 | 1 | 0 | 0 | 1 | 1 | 1 | 2 | 1 | 1 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1.0 | 20 | 9 |
| 3 | 1 | 1 | 0 | 0 | 1 | 1 | 1 | 6 | 3 | 2 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 2 | 0.5 | 36 | 8 |
| 4 | 1 | 1 | 0 | 0 | 1 | 1 | 2 | 1 | 2 | 3 | ... | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1.0 | 36 | 13 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 1762237 | 9 | 7 | 1 | 0 | 2 | 35 | 2 | 1 | 1 | 1 | ... | 0 | 0 | 0 | 0 | 0 | 4 | 1 | 1.0 | 52 | 5 |
| 1762238 | 9 | 7 | 1 | 0 | 2 | 36 | 2 | 1 | 1 | 1 | ... | 0 | 0 | 0 | 0 | 0 | 1 | 3 | 1.0 | 27 | 6 |
| 1762239 | 9 | 7 | 1 | 0 | 2 | 36 | 2 | 1 | 1 | 1 | ... | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 0.5 | 39 | 10 |
| 1762240 | 9 | 7 | 1 | 0 | 2 | 4 | 2 | 1 | 1 | 1 | ... | 0 | 0 | 0 | 0 | 0 | 1 | 3 | 1.0 | 41 | 10 |
| 1762241 | 9 | 7 | 1 | 0 | 5 | 24 | 1 | 1 | 1 | 1 | ... | 0 | 0 | 0 | 0 | 0 | 1 | 5 | 1.0 | 44 | 10 |
1762242 rows × 34 columns
data.isna().sum()
C_MNTH 0 C_WDAY 0 C_HOUR 0 C_SEV 0 C_VEHS 0 C_CONF 0 C_RCFG 0 C_WTHR 0 C_RSUR 0 C_RALN 0 C_TRAF 0 V_TYPE_1 0 V_TYPE_10 0 V_TYPE_11 0 V_TYPE_14 0 V_TYPE_16 0 V_TYPE_17 0 V_TYPE_18 0 V_TYPE_19 0 V_TYPE_20 0 V_TYPE_21 0 V_TYPE_22 0 V_TYPE_23 0 V_TYPE_24 0 V_TYPE_5 0 V_TYPE_6 0 V_TYPE_7 0 V_TYPE_8 0 V_TYPE_9 0 P_SEX_F 0 P_SEX_M 0 P_SAFE 0 P_AGE 0 A_VAGE 0 dtype: int64
target = ['C_SEV']
v_type_col = [col for col in data.columns if col[0:6]=="V_TYPE"]
p_sex_col = [col for col in data.columns if col[0:5]=="P_SEX"]
int_features = ['C_VEHS', 'P_AGE', 'A_VAGE'] + v_type_col + p_sex_col
float_features = ['P_SAFE']
numeric_features = int_features + float_features
categories = [col for col in data.columns if col not in target+numeric_features]
dtypes = {}
for feature in int_features:
dtypes[feature] = 'uint8'
for feature in float_features:
dtypes[feature] = 'float'
for feature in categories+target:
dtypes[feature] = 'uint8'
data = data.astype(dtypes)
data[categories] = data[categories].astype('category')
data.dtypes.to_dict()
{'C_MNTH': CategoricalDtype(categories=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12], ordered=False),
'C_WDAY': CategoricalDtype(categories=[1, 2, 3, 4, 5, 6, 7], ordered=False),
'C_HOUR': CategoricalDtype(categories=[0, 1, 2, 3], ordered=False),
'C_SEV': dtype('uint8'),
'C_VEHS': dtype('uint8'),
'C_CONF': CategoricalDtype(categories=[1, 2, 3, 4, 5, 6, 21, 22, 23, 24, 25, 31, 32, 33, 34, 35, 36,
41, 42],
, ordered=False),
'C_RCFG': CategoricalDtype(categories=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 13], ordered=False),
'C_WTHR': CategoricalDtype(categories=[1, 2, 3, 4, 5, 6, 7, 8], ordered=False),
'C_RSUR': CategoricalDtype(categories=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10], ordered=False),
'C_RALN': CategoricalDtype(categories=[1, 2, 3, 4, 5, 6, 7], ordered=False),
'C_TRAF': CategoricalDtype(categories=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 15, 16, 17, 18,
19],
, ordered=False),
'V_TYPE_1': dtype('uint8'),
'V_TYPE_10': dtype('uint8'),
'V_TYPE_11': dtype('uint8'),
'V_TYPE_14': dtype('uint8'),
'V_TYPE_16': dtype('uint8'),
'V_TYPE_17': dtype('uint8'),
'V_TYPE_18': dtype('uint8'),
'V_TYPE_19': dtype('uint8'),
'V_TYPE_20': dtype('uint8'),
'V_TYPE_21': dtype('uint8'),
'V_TYPE_22': dtype('uint8'),
'V_TYPE_23': dtype('uint8'),
'V_TYPE_24': dtype('uint8'),
'V_TYPE_5': dtype('uint8'),
'V_TYPE_6': dtype('uint8'),
'V_TYPE_7': dtype('uint8'),
'V_TYPE_8': dtype('uint8'),
'V_TYPE_9': dtype('uint8'),
'P_SEX_F': dtype('uint8'),
'P_SEX_M': dtype('uint8'),
'P_SAFE': dtype('float64'),
'P_AGE': dtype('uint8'),
'A_VAGE': dtype('uint8')}
En primer lugar, se comprueba el número de ceros y unos de la variable objetivo para decidir si se debe hacer un oversampling o no. Mientras que los accidentes sin fatalidades exceden el millón ochocientos mil, los que provocan muertes ni siquiera llegan a cuarenta mil. Por lo tanto, se procede a la realización de un oversampling que eleve la muestra de unos.
Para ello, se divide la muestra primeramente en train y test para que la división mantenga las características originales del dataset.
data_major = data[data.C_SEV == 0]
data_minor = data[data.C_SEV == 1]
print(data_major)
print(data_minor)
C_MNTH C_WDAY C_HOUR C_SEV C_VEHS C_CONF C_RCFG C_WTHR C_RSUR \
0 1 1 0 0 1 1 1 1 4
1 1 1 0 0 1 1 1 1 5
2 1 1 0 0 1 1 1 2 1
3 1 1 0 0 1 1 1 6 3
4 1 1 0 0 1 1 2 1 2
... ... ... ... ... ... ... ... ... ...
1762237 9 7 1 0 2 35 2 1 1
1762238 9 7 1 0 2 36 2 1 1
1762239 9 7 1 0 2 36 2 1 1
1762240 9 7 1 0 2 4 2 1 1
1762241 9 7 1 0 5 24 1 1 1
C_RALN ... V_TYPE_5 V_TYPE_6 V_TYPE_7 V_TYPE_8 V_TYPE_9 P_SEX_F \
0 1 ... 0 1 0 0 0 1
1 1 ... 0 0 0 0 0 0
2 1 ... 0 0 0 0 0 0
3 2 ... 0 0 0 0 0 0
4 3 ... 0 0 0 0 0 0
... ... ... ... ... ... ... ... ...
1762237 1 ... 0 0 0 0 0 4
1762238 1 ... 0 0 0 0 0 1
1762239 1 ... 0 0 0 0 0 1
1762240 1 ... 0 0 0 0 0 1
1762241 1 ... 0 0 0 0 0 1
P_SEX_M P_SAFE P_AGE A_VAGE
0 0 1.0 17 8
1 1 0.0 30 18
2 1 1.0 20 9
3 2 0.5 36 8
4 1 1.0 36 13
... ... ... ... ...
1762237 1 1.0 52 5
1762238 3 1.0 27 6
1762239 1 0.5 39 10
1762240 3 1.0 41 10
1762241 5 1.0 44 10
[1726427 rows x 34 columns]
C_MNTH C_WDAY C_HOUR C_SEV C_VEHS C_CONF C_RCFG C_WTHR C_RSUR \
116 1 1 1 1 2 21 2 4 3
117 1 1 1 1 2 31 1 4 3
118 1 1 1 1 2 32 1 4 4
209 1 1 1 1 2 31 1 1 5
320 1 1 2 1 1 1 2 1 1
... ... ... ... ... ... ... ... ... ...
1762083 9 7 3 1 2 21 1 1 1
1762084 9 7 3 1 3 32 1 1 1
1762109 9 7 3 1 2 35 2 1 1
1762147 9 7 0 1 1 4 2 1 1
1762171 9 7 0 1 2 21 1 1 1
C_RALN ... V_TYPE_5 V_TYPE_6 V_TYPE_7 V_TYPE_8 V_TYPE_9 P_SEX_F \
116 2 ... 0 0 1 0 0 1
117 3 ... 0 0 0 0 0 2
118 4 ... 0 0 0 1 0 0
209 4 ... 0 0 0 1 0 1
320 1 ... 0 0 0 0 0 0
... ... ... ... ... ... ... ... ...
1762083 1 ... 0 0 3 0 0 1
1762084 1 ... 0 0 0 0 0 3
1762109 1 ... 0 0 0 0 0 1
1762147 4 ... 0 0 0 0 0 0
1762171 1 ... 0 0 0 1 0 1
P_SEX_M P_SAFE P_AGE A_VAGE
116 2 1.000000 58 5
117 1 1.000000 56 4
118 3 0.666667 69 12
209 3 1.000000 27 11
320 2 1.000000 49 6
... ... ... ... ...
1762083 2 1.000000 63 21
1762084 3 0.833333 42 2
1762109 2 0.666667 49 8
1762147 1 1.000000 25 35
1762171 2 0.666667 23 15
[35815 rows x 34 columns]
pd_plot_target = data['C_SEV']\
.value_counts(normalize=True)\
.mul(100).rename('percent').reset_index()
pd_plot_target_conteo = data['C_SEV'].value_counts().reset_index()
pd_plot_target_pc = pd.merge(pd_plot_target, pd_plot_target_conteo, on=['index'], how='inner')
fig = px.histogram(pd_plot_target_pc, x="index", y=['percent'])
fig.show()
X_arr = np.array(data.drop('C_SEV',axis=1))
y_arr = np.array(data['C_SEV'])
X_train, X_test, y_train, y_test = train_test_split(X_arr, y_arr, test_size=0.2,
random_state=12345)
X_train = pd.DataFrame(X_train, columns=[col for col in data.columns if col not in target])
dtypes = {}
for feature in int_features:
dtypes[feature] = 'uint8'
for feature in float_features:
dtypes[feature] = 'float'
for feature in categories:
dtypes[feature] = 'uint8'
X_train = X_train.astype(dtypes)
X_test = pd.DataFrame(X_test, columns=[col for col in data.columns if col not in target])
dtypes = {}
for feature in int_features:
dtypes[feature] = 'uint8'
for feature in float_features:
dtypes[feature] = 'float'
for feature in categories:
dtypes[feature] = 'uint8'
X_test = X_test.astype(dtypes)
y_test = pd.DataFrame(y_test, columns=[col for col in data.columns if col in target])
y_test = y_test.astype('uint8')
y_train = pd.DataFrame(y_train, columns=[col for col in data.columns if col in target])
y_train = y_train.astype('uint8')
data_train = pd.concat([X_train, y_train],axis=1)
data_test = pd.concat([X_test, y_test],axis=1)
print('== Train\n', data_train['C_SEV'].value_counts(normalize=True))
print('== Test\n', data_test['C_SEV'].value_counts(normalize=True))
== Train 0 0.979627 1 0.020373 Name: C_SEV, dtype: float64 == Test 0 0.979875 1 0.020125 Name: C_SEV, dtype: float64
Con la función SMOTE se realiza un oversampling automático que necesita de una semilla (random_state) para obtener los valores requeridos y de un sampling_strategy que determine a grandes rasgos el porcentaje al que queremos llegar con la ampliación del muestreo. Cabe destacar que este procedimiento se realiza únicamente sobre el training, dejando el test aparte con las propiedades orginales del dataset
sm = SMOTE(random_state=12345, n_jobs=-1, sampling_strategy=0.50)
X_train_oversampled, y_train_oversampled = sm.fit_resample(X_train, y_train)
X_train_oversampled = pd.DataFrame(X_train_oversampled, columns=X_train.columns)
print("Distribution before resampling {}".format(Counter(Y_train)))
print("Distribution after resampling {}".format(Counter(Y_train_oversampled)))
Distribution before resampling Counter({0: 1381071, 1: 28722})
Distribution after resampling Counter({0: 1381071, 1: 690535})
pd_plot_target_pc = Y_train_oversampled\
.value_counts(normalize=True)\
.mul(100).rename('percent').reset_index()
pd_plot_target_conteo = Y_train_oversampled.value_counts().reset_index()
pd_plot_target = pd.merge(pd_plot_target_pc,
pd_plot_target_conteo, on='index', how='inner')
fig = px.histogram(pd_plot_target, x="index", y=['percent'])
fig.show()
En primer lugar, se realizará un modelado del dataset con la configuración por defecto, para tener una primera referencia de la performance de los dstintos modelos. Para ello, hacemos una lista de los modelos que queremos estudiar: GLM o modelo lineal generalizado (Tweedie Regressor en nuestro caso), Random Forest, XGB y LightGBM. Posteriormente, con un bucle y el empleo del Pipeline podemos estandarizar los valores de cada variable para cada clasificador y entrenar así todos los modelos. Una vez ejecutados, se comprueba que los que mejor score tienen son el Random Forest y el LightGBM, ambos con un 98%.
preprocessor = ColumnTransformer(
transformers=[
('scaler', StandardScaler(), numeric_features),
('te', ce.TargetEncoder(), categories)], n_jobs=-1)
classifiers = [
linear_model.LogisticRegression(penalty='none', n_jobs=-1),
RandomForestClassifier(n_jobs=-1),
XGBClassifier(objective="binary:logistic", n_jobs=-1),
LGBMClassifier(objective="binary" , n_jobs=-1)
]
for classifier in classifiers:
pipe = Pipeline(steps=[('preprocessor', preprocessor),
('classifier', classifier)])
pipe.fit(X_train_oversampled, y_train_oversampled)
print(classifier)
print(evaluate_model(y_test, pipe.predict(X_test)))
LogisticRegression(n_jobs=-1, penalty='none')
Accuracy: 0.89003231673235
Recall: 0.3335683067813337
Precision: 0.06500535758441629
F-Score: 0.10881
Confusion matrix:
[[311325 34031]
[ 4727 2366]]
None
RandomForestClassifier(n_jobs=-1)
Accuracy: 0.9747821670653058
Recall: 0.09248554913294797
Precision: 0.2111361441905375
F-Score: 0.12863
Confusion matrix:
[[342905 2451]
[ 6437 656]]
None
[15:06:30] WARNING: ..\src\learner.cc:1115: Starting in XGBoost 1.3.0, the default evaluation metric used with the objective 'binary:logistic' was changed from 'error' to 'logloss'. Explicitly set eval_metric if you'd like to restore the old behavior.
XGBClassifier(base_score=0.5, booster='gbtree', colsample_bylevel=1,
colsample_bynode=1, colsample_bytree=1, enable_categorical=False,
gamma=0, gpu_id=-1, importance_type=None,
interaction_constraints='', learning_rate=0.300000012,
max_delta_step=0, max_depth=6, min_child_weight=1, missing=nan,
monotone_constraints='()', n_estimators=100, n_jobs=-1,
num_parallel_tree=1, predictor='auto', random_state=0,
reg_alpha=0, reg_lambda=1, scale_pos_weight=1, subsample=1,
tree_method='exact', validate_parameters=1, verbosity=None)
Accuracy: 0.9499303445321167
Recall: 0.22924009586916677
Precision: 0.11777488048674489
F-Score: 0.15561
Confusion matrix:
[[333176 12180]
[ 5467 1626]]
None
LGBMClassifier(objective='binary')
Accuracy: 0.9497941546152777
Recall: 0.23628929930917805
Precision: 0.12010892933925756
F-Score: 0.15926
Confusion matrix:
[[333078 12278]
[ 5417 1676]]
None
Como cualquier modelización mediante regresión lineal, es conveniente escalar las variables antes de estimar el modelo, ya que el algoritmo es sensible a la magnitud de los datos, y puede influir negativamente en la estimación de los parámetros. Por ello, este será el primer paso a realizar.
pipe_lasso = Pipeline(steps=[('preprocessor', preprocessor),
('lasso_model', LogisticRegression(C=1, penalty='l1',
solver='liblinear', n_jobs=-1))])
pipe_lasso.fit(X_train_oversampled, y_train_oversampled)
coefficients = pipe_lasso.named_steps['lasso_model'].coef_
importance = np.abs(coefficients)
preprocessor = pipe_lasso.named_steps["preprocessor"]
features = get_feature_names(preprocessor)
df_coeficientes_lasso = pd.DataFrame({'predictor':coefficients.flatten(),
'coef': features}) #index=
fig, ax = plt.subplots(figsize=(16, 3.84))
ax.stem(df_coeficientes_lasso.coef, df_coeficientes_lasso.predictor, markerfmt=' ')
plt.xticks(rotation=90, ha='right', size=10)
ax.set_xlabel('variable')
ax.set_ylabel('coeficientes')
ax.set_title('Coeficientes del modelo lasso');
plt.show()
Aquellas variables que obtienen una importancia mayor tras la penalización Lasso son la V_TYPE_1 (coches), la V_TYPE_17 (bicicletas), la P_SEX_M (hombres) y la P_SAFE (porcentaje de personas involucradas en el accidente que utilizaban algún tipo de dispositivo de seguridad)
aux=pd.DataFrame(importance.tolist(), columns=features)
aux = aux.T
aux.columns = ['importance']
selected_feat_lasso = aux[aux['importance'] >= 0.05] #threshold
print('total features: {}'.format((len(features))))
print('selected features: {}'.format(len(selected_feat_lasso)))
total features: 33 selected features: 28
A continuación se va a llevar a cabo un nuevo entrenamiento de los clasficiadores utilizados en la primera modelización, esta vez considerando únicamente aquellas variables significativas tras la regulariazación Lasso. Por lo tanto, el primer paso será modificar el conjunto de datos de X_train y X_test, para que contengan únicamente las variables importantes.
selected_feat_lasso = selected_feat_lasso.index.values.tolist()
features_lasso = [i.split('__')[1] for i in selected_feat_lasso]
features_lasso
['C_VEHS', 'P_AGE', 'A_VAGE', 'V_TYPE_1', 'V_TYPE_10', 'V_TYPE_11', 'V_TYPE_14', 'V_TYPE_16', 'V_TYPE_17', 'V_TYPE_19', 'V_TYPE_20', 'V_TYPE_22', 'V_TYPE_23', 'V_TYPE_24', 'V_TYPE_5', 'V_TYPE_6', 'V_TYPE_7', 'V_TYPE_8', 'V_TYPE_9', 'P_SEX_F', 'P_SEX_M', 'P_SAFE', 'C_HOUR', 'C_RCFG', 'C_WTHR', 'C_RSUR', 'C_RALN', 'C_TRAF']
X_train_short = X_train_oversampled[features_lasso]
X_test_short = X_test[features_lasso]
numeric_features_short = [col for col in features_lasso if col in numeric_features]
categories_short = [col for col in features_lasso if col in categories]
preprocessor = ColumnTransformer(
transformers=[
('scaler', StandardScaler(), numeric_features_short),
('ohe', ce.TargetEncoder(), categories_short)], n_jobs=-1)
for classifier in classifiers:
pipe = Pipeline(steps=[('preprocessor', preprocessor),
('classifier', classifier)])
pipe.fit(X_train_short, y_train_oversampled)
print(classifier)
print(evaluate_model(y_test, pipe.predict(X_test_short)))
LogisticRegression(n_jobs=-1, penalty='none')
Accuracy: 0.8908040595944378
Recall: 0.3299027209925278
Precision: 0.06486846117594877
F-Score: 0.10842
Confusion matrix:
[[311623 33733]
[ 4753 2340]]
None
RandomForestClassifier(n_jobs=-1)
Accuracy: 0.9598778830412343
Recall: 0.0843084731425349
Precision: 0.0725376031052887
F-Score: 0.07798
Confusion matrix:
[[337710 7646]
[ 6495 598]]
None
[20:59:00] WARNING: ..\src\learner.cc:1115: Starting in XGBoost 1.3.0, the default evaluation metric used with the objective 'binary:logistic' was changed from 'error' to 'logloss'. Explicitly set eval_metric if you'd like to restore the old behavior.
XGBClassifier(base_score=0.5, booster='gbtree', colsample_bylevel=1,
colsample_bynode=1, colsample_bytree=1, enable_categorical=False,
gamma=0, gpu_id=-1, importance_type=None,
interaction_constraints='', learning_rate=0.300000012,
max_delta_step=0, max_depth=6, min_child_weight=1, missing=nan,
monotone_constraints='()', n_estimators=100, n_jobs=-1,
num_parallel_tree=1, predictor='auto', random_state=0,
reg_alpha=0, reg_lambda=1, scale_pos_weight=1, subsample=1,
tree_method='exact', validate_parameters=1, verbosity=None)
Accuracy: 0.9505715720572339
Recall: 0.15776117298745243
Precision: 0.08904981696641731
F-Score: 0.11384
Confusion matrix:
[[333909 11447]
[ 5974 1119]]
None
LGBMClassifier(objective='binary')
Accuracy: 0.9523136680767997
Recall: 0.16692513745946708
Precision: 0.09799702036086741
F-Score: 0.12349
Confusion matrix:
[[334458 10898]
[ 5909 1184]]
None
Por lo tanto, aquel clasificador que se presenta como el que mejor estimación del modelo proporciona es el LightGBM con unas métricas superiores al resto de clasificadores destacando principalmente en su Recall y F-Score respecto a los demás. Además, observando la matrices de confusión, es aquel que es capaz de identificar un mayor número de valor true positive. Por lo tanto, será aquel utilizado para optimizar a través de la selección de hiperparámetros mediante CrossValidation. El XGBoost, por su parte, tiene una performance similar, siendo ligeramente inferior. EL LinearModel LogisticRegression, es capaz de detectar un gran número de valores TP, pero sin embargo comete un mayor error que perjudica su performance, como muestran sus principales métricas. Por último, cabe destacar que el clasificador RandomForest se presenta como el que peor es capaz de estimar el modelo. Además, se observa que el rendimiento de todos los modelos es inferior tras llevar a cabo la reducción de variables mediante la estimación de un modelo lineal Lasso con penalización L1, por lo que se considerarán todas las variables en la optimización por hiperparámetros. Esto tiene sentido, ya que un mayor número de variables otorga una mayor flexibilidad a los modelos basados en árboles y les permite tener un mejor comportamiento.
data_train = pd.concat([y_train_oversampled, X_train_oversampled], axis=1)
data_test = pd.concat([y_test, X_test], axis=1)
# Guardar el dataset transformado
data_train.to_csv("data/pd_data_train.csv", index=False)
data_test.to_csv("data/pd_data_test.csv", index=False)
Para la optimización de hiperparámetros del modelo seleccionado anteriormente, se va a utilizar la técnica de validación cruzada. Se ha decidido implementar dos fases en el proceso. Una primera fase, se utilizará el RandomizedSearchCV para acotar el espacio de hiperprámetros. Con aquellos parámetros seleccionados, se realizará una nueva validación cruzada utilizando GridSearchCV sobre valores cercanos a los hiperparámetros devueltos en la primera fase, de manera que se realice una búsqueda más exhaustiva hasta llegar a la configuración óptima para el problema.
preprocessor = ColumnTransformer(
transformers=[
('scaler', StandardScaler(), numeric_features),
('te', ce.TargetEncoder(), categories)], n_jobs=-1)
pipe_lgbmc_ini = Pipeline(steps=[('preprocessor', preprocessor),
('classifier', LGBMClassifier())])
### Parameter Tunning Optimization
param_grid = {
'classifier__objective': ["binary"],
'classifier__boosting': ["gbdt", "rf", "dart", "goss"],
'classifier__n_estimators':[100,200],
'classifier__num_leaves' : [31,64,128,200],
'classifier__max_depth' :[-1,5,7,8,10],
'classifier__learning_rate': [0.08, 0.1, 0.2, 0.5, 0.8, 1]}
# this will train 100 models over 5 folds of cross validation (500 models total)
CV = RandomizedSearchCV(pipe_lgbmc_ini, param_grid, cv=10, n_jobs=-1, random_state=12345) #n_iter=100, random_state=12345,
CV.fit(X_train_oversampled, y_train_oversampled)
print(CV.best_params_)
print(CV.best_score_)
[LightGBM] [Warning] boosting is set=dart, boosting_type=gbdt will be ignored. Current value: boosting=dart
{'classifier__objective': 'binary', 'classifier__num_leaves': 200, 'classifier__n_estimators': 200, 'classifier__max_depth': -1, 'classifier__learning_rate': 1, 'classifier__boosting': 'dart'}
0.9297989318899443
preprocessor = ColumnTransformer(
transformers=[
('scaler', StandardScaler(), numeric_features),
('te', ce.TargetEncoder(), categories)], n_jobs=-1)
pipe_lgbmc = Pipeline(steps=[('preprocessor', preprocessor),
('classifier', LGBMClassifier())])
### Parameter Tunning Optimization
param_grid = {
'classifier__objective': ["binary"],
'classifier__boosting': ["gbdt", "dart"],
'classifier__n_estimators':[200, 250],
'classifier__num_leaves' : [180, 200, 250],
'classifier__max_depth' :[-1],
'classifier__learning_rate': [0.1, 1]}
# this will train 100 models over 5 folds of cross validation (500 models total)
CV = GridSearchCV(pipe_lgbmc, param_grid, cv=10, n_jobs=-1) #n_iter=100, random_state=12345,
CV.fit(X_train_oversampled, y_train_oversampled)
print(CV.best_params_)
print(CV.best_score_)
{'classifier__boosting': 'gbdt', 'classifier__learning_rate': 0.1, 'classifier__max_depth': -1, 'classifier__n_estimators': 250, 'classifier__num_leaves': 250, 'classifier__objective': 'binary'}
0.9383690900647554
predictions = CV.predict(X_test)
evaluate_model(y_test, predictions)
Accuracy: 0.9598126253727489 Recall: 0.19061046101790496 Precision: 0.13831202046035806 F-Score: 0.16030 Confusion matrix: [[336933 8423] [ 5741 1352]]
titles_options = [("Confusion matrix, without normalization", None),
("Normalized confusion matrix", 'true')]
for title, normalize in titles_options:
disp = plot_confusion_matrix(CV, X_test, y_test,
# display_labels=ytest,
cmap=plt.cm.Blues,
normalize=normalize)
disp.ax_.set_title(title)
print(title)
print(disp.confusion_matrix)
plt.show()
Confusion matrix, without normalization [[336933 8423] [ 5741 1352]] Normalized confusion matrix [[0.97561067 0.02438933] [0.80938954 0.19061046]]
prob_predictions = CV.predict_proba(X_test)
yhat = prob_predictions[:, 1]
fpr, tpr, thresholds = roc_curve(y_test, yhat)
# calculate the g-mean for each threshold
gmeans = np.sqrt(tpr * (1-fpr))
# locate the index of the largest g-mean
ix = np.argmax(gmeans)
print('Best Threshold=%f, G-Mean=%.3f' % (thresholds[ix], gmeans[ix]))
Best Threshold=0.075053, G-Mean=0.721
Se obtiene que el threshold que optimiza el modelo es igual a 0.075053
# plot the roc curve for the model
plt.plot([0,1], [0,1], linestyle='--', label='No Skill')
plt.plot(fpr, tpr, marker='.', label='LigtGBM')
plt.scatter(fpr[ix], tpr[ix], s=100, marker='o', color='black', label='Best')
# axis labels
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.legend()
# show the plot
plt.show()
La Curva de ROC cuenta con varios parámetros que la determinan. El primero, la g-mean (72.1%) que mide la bondad de acierto tanto de 1 como de 0 que tiene el modelo y se corresponde con el área bajo la curva (aciertos). Todos los puntos que se sitúan por encima de la curva se "escapan al modelo", ya que no se identifican correctamente. A su vez, el gráfico muestra el threshold óptimo para nuestro problema, que permite maximizar la tasa de acierto de T minimizando la tasa de fallo de FP. En este caso, ese valor sería igual a 0.075053, de forma que en las preddicciones en que se obtenga una probabilidad de ocurrencia de fatalidad superior a ese número, pasarían a ser considerados como unos y todos los de debajo como ceros. En un principio el threshold considerado para llevar a cabo dicha clasificación es de 0.5. No obstante, como nuestro modelo detecta peor los unos, tiene sentido que el umbral de clasificación baje para que aumente la detección, en detrimento de un mejor porcentaje de acierto en los ceros.
# calculate pr-curve
precision, recall, thresholds = precision_recall_curve(y_test, yhat)
# convert to f score
fscore = (2 * precision * recall) / (precision + recall)
# locate the index of the largest f score
ix = np.argmax(fscore)
print('Best Threshold=%f, F-Score=%.3f' % (thresholds[ix], fscore[ix]))
# plot the roc curve for the model
no_skill = len(y_test[y_test==1]) / len(y_test)
plt.plot([0,1], [no_skill,no_skill], linestyle='--', label='No Skill')
plt.plot(recall, precision, marker='.', label='LightGBM')
plt.scatter(recall[ix], precision[ix], s=100, marker='o', color='black', label='Best')
# axis labels
plt.xlabel('Recall')
plt.ylabel('Precision')
plt.legend()
# show the plot
plt.show()
Best Threshold=0.975579, F-Score=nan
skplt.metrics.plot_cumulative_gain(y_test, prob_predictions)
plt.show()
Como se puede observar en el gráfico de la función acumulativa de ganancia como consecuenia de que el dataset está muy desbalanceado, el punto de ganancia máximo se alcanza alrededor del percentil 20%, es decir, que con el 20% de accidentes con mayor probabilidad de fatalidad, se obtiene el 60% del total de accidentes con fatalidad.
skplt.metrics.plot_lift_curve(y_test, prob_predictions)
plt.show()
Como se puede observar en el gráfico de la curva lift como consecuenia de que el dataset está muy desbalanceado, la curva es muy plana, siendo aquellos percentiles con una probabilidad mayor los que obtienen un mayor lift, que es el ratio que compara el porcentaje de aciertos de TP en ese punto respecto el porcentaje de aciertos en el total de la muestra.
Acontinuación, se realizarán varias pruebas en las que se modificará el valor del threshold para detectar el modelo que se ajuste mejor a los requerimientos del problema. Comenzamos por el Threshold óptimo calculado anteriormente:
y_pred_best = (CV.predict_proba(X_test)[:,1] >= 0.075053).astype(bool)
y_pred_best
array([ True, False, True, ..., False, False, True])
evaluate_model(y_test, y_pred_best)
Accuracy: 0.7293934725307775 Recall: 0.7125334837163401 Precision: 0.051367008842362026 F-Score: 0.09583 Confusion matrix: [[252020 93336] [ 2039 5054]]
cm_gbt = confusion_matrix(y_test, y_pred_best, normalize='true')
print('')
plt.figure(figsize=(6,6))
sns.heatmap(cm_gbt, annot=True, fmt='.5f', linewidths=.5, square = True, cmap = 'Blues_r');
plt.ylabel('Actual label')
plt.show();
y_pred_best_2 = (CV.predict_proba(X_test)[:,1] >= 0.15).astype(bool)
y_pred_best_2
array([ True, False, True, ..., False, False, False])
evaluate_model(y_test, y_pred_best_2)
Accuracy: 0.8457734310495987 Recall: 0.5384181587480614 Precision: 0.06956030745692325 F-Score: 0.12320 Confusion matrix: [[294273 51083] [ 3274 3819]]
cm_gbt = confusion_matrix(y_test, y_pred_best_2, normalize='true')
print('')
plt.figure(figsize=(6,6))
sns.heatmap(cm_gbt, annot=True, fmt='.5f', linewidths=.5, square = True, cmap = 'Blues_r');
plt.ylabel('Actual label')
plt.show();
y_pred_best_3 = (CV.predict_proba(X_test)[:,1] >= 0.12).astype(bool)
y_pred_best_3
array([ True, False, True, ..., False, False, True])
evaluate_model(y_test, y_pred_best_3)
Accuracy: 0.7733033715516288 Recall: 0.6297758353306077 Precision: 0.054648886713971126 F-Score: 0.10057 Confusion matrix: [[268083 77273] [ 2626 4467]]
cm_gbt = confusion_matrix(y_test, y_pred_best_3, normalize='true')
print('')
plt.figure(figsize=(6,6))
sns.heatmap(cm_gbt, annot=True, fmt='.5f', linewidths=.5, square = True, cmap = 'Blues_r');
plt.ylabel('Actual label')
plt.show();
| Ratios | Sin Selecc.Var | Selecc. var. | CV | TH optimo |
|---|---|---|---|---|
| Accuracy | 0.94979 | 0.95231 | 0.95981 | 0.72939 |
| Recall | 0.23629 | 0.16693 | 0.19061 | 0.71253 |
| Precision | 0.12011 | 0.09799 | 0.13831 | 0.05137 |
| F-Score | 0.15926 | 0.12349 | 0.16030 | 0.09583 |
En la tabla anterior se detalla la evolución de las principales métricas de evaluación de los modelos estimados. El primer modelo, lanzado con la configuración de hiperparámetros por defecto, muestra un rendimiento correcto, con un accuracy del 95% y un F1-Score del 16% como principales métricas. Tras la reducción de variables con la estimación de un modelo lineal Lasso con penalización L1, se observa que el performance del modelo se reduce significativamente, disminuyendo todas las métricas a excepción del Accuracy, por lo que en la optimización posterior de hiperparámetros se han considerado todas las variables. La estimación del modelo LightGBM tras la optimización de sus parámetros mediante validación cruzada, mejora el rendimiento del modelo, pero no lo hace de manera significativa. Esto, puede en parte producirse, por algo de overfitting ya que tanto el número de árboles estimados como el número de hojas seleccionado era el mayor posible. Tras realizar varias comprobaciones variando el threshold del último modelo estimado, se ha determinado que el threshold óptimo es aquel que permite modelizar mejor la serie de datos para el objetivo planteado, ya que aunque se reduce tanto el Accuracy como el Precission como el F1-Score, es aquel que devuelve una matriz de confusión más balanceada en cuanto al número de TP y TN, y como consecuencia aquel que permite maximizar el Recall. El Recall, para el caso de estudio se considera una métrica clave, ya que no predecir correctamente una fatalidad y que realmente exista puede ser más perjudicial para la aseguradora que el caso de predecir fatalidad y que finalmente no ocurra. Por lo tanto, el modelo seleccionado, con un threshold de 0.075053, predecirá fatalidad cuando la pred_proba de 1 sea mayor o igual al 7,5%. Que el threshold óptimo sea tan pequeño, se produce como consecuencia de las características del dataset que está muy desbalanceado. Por último destacar como la matriz de confusión del modelo arroja unos resultados de acierto del 73% en el caso de los accidentes sin fatalidad y del 71% de los accidentes que suponen al menos una fatalidad.
La interpretabilidad es el grado en el que una persona puede predecir el resultado de un modelo. Cuanta mayor interpretabilidad tenga un modelo, mas fácil será comprenderlo y predecirlo. Es una característica necesaria para paliar las situaciones en las que la información está incompleta, ya que para ciertos proyectos o problemas la predicción por sí sola no es suficiente, sino que necesita de una explicación de cómo se ha llegado a ella: interpretabilidad del modelo. A su vez, la interpretación de una predicción errónea ayuda a comprender la causa del error. Los principales rasgos que todo modelo con una buena interpretabilidad ha de cumplir son los siguientes:
# load JS visualization code to notebook
shap.initjs()
explainer = shap.TreeExplainer(CV.best_estimator_.named_steps['classifier'])
preprocessor = CV.best_estimator_.named_steps["preprocessor"] #pipe_lgbmc.
test_data = pd.DataFrame(CV.best_estimator_.named_steps['preprocessor'].transform(X_test),
columns=X_train.columns.tolist()) #fit_ (, y_test) .sparse.from_spmatrix
test_data.head()
| C_MNTH | C_WDAY | C_HOUR | C_VEHS | C_CONF | C_RCFG | C_WTHR | C_RSUR | C_RALN | C_TRAF | ... | V_TYPE_5 | V_TYPE_6 | V_TYPE_7 | V_TYPE_8 | V_TYPE_9 | P_SEX_F | P_SEX_M | P_SAFE | P_AGE | A_VAGE | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | -0.968343 | -1.315754 | -0.559305 | -0.203803 | -0.007751 | -0.050718 | -0.209312 | -0.054875 | -0.19696 | -0.026587 | ... | 0.528139 | 9.0 | 1.0 | 2.0 | 6.0 | 1.0 | 2.0 | 1.0 | 2.0 | 18.0 |
| 1 | -0.968343 | 0.720556 | -0.366535 | -0.975164 | -0.007751 | -0.050718 | -0.209312 | -0.054875 | -0.19696 | -0.026587 | ... | 0.528139 | 11.0 | 3.0 | 1.0 | 4.0 | 1.0 | 2.0 | 5.0 | 2.0 | 18.0 |
| 2 | 0.423383 | -0.772738 | -0.270150 | 0.181878 | -0.007751 | -0.050718 | -0.209312 | -0.054875 | -0.19696 | -0.026587 | ... | 0.528139 | 1.0 | 5.0 | 1.0 | 31.0 | 1.0 | 1.0 | 3.0 | 3.0 | 18.0 |
| 3 | -0.968343 | 0.177540 | -0.270150 | -0.589483 | -0.007751 | -0.050718 | -0.209312 | -0.054875 | -0.19696 | -0.026587 | ... | 0.528139 | 4.0 | 4.0 | 0.0 | 2.0 | 2.0 | 1.0 | 1.0 | 2.0 | 3.0 |
| 4 | 0.423383 | 0.041786 | -0.366535 | -0.203803 | -0.007751 | -0.050718 | -0.209312 | -0.054875 | -0.19696 | -0.026587 | ... | 0.528139 | 11.0 | 2.0 | 2.0 | 35.0 | 2.0 | 1.0 | 1.0 | 1.0 | 1.0 |
5 rows × 33 columns
shap_values = explainer.shap_values(test_data)
len(shap_values[0][0])
33
# visualize the first prediction's explanation (use matplotlib=True to avoid Javascript)
shap.force_plot(explainer.expected_value[0], shap_values[0][0], test_data.iloc[0,:], matplotlib=True)
shap.force_plot(explainer.expected_value[0], shap_values[0][:100,:], features=test_data.iloc[:100, :])
shap.summary_plot(shap_values, X_test)
Como se observa en el gráfico, ninguna de las variables demuestra un elevado grado discriminatorio sobre la variable objetivo.
shap.summary_plot(shap_values[0], features=X_test, max_display=10)
En el gráfico se representan las 10 variables con mayor importancia. Por lo tanto las principales variables que afectan al modelo tienen una contribución marginal media positiva y tienden al valor 0. A excepción de la variable P_SAFE, que discrimina en mayor medida hacia el valor 1. En un punto intermedio se encuentran el C_VEHS y la A_VAGE.
Conclusiones de la interpretabilidad: